Efficient mixed-force first-principles molecular dynamics 
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We present an efficient method to mix well converged ab initio forces with simpler and faster ones 
in molecular dynamics. While the cheap forces are evaluated every time step, the converged ones 
correct the trajectory only every n time steps. For convenience, both types of forces are calculated 
with the same basic scheme, using density functional theory, norm conserving pseudopotentials, and 
a basis set of numerical atomic orbitals. The cheap forces are evaluated with a short-range minimal 
basis set and the non-selfconsistent Harris functional. Since these evaluations are hundreds of times 
faster than those of the converged forces, they add a neglegible cost, and the boost in computational 
efficiency is approximately a factor n. Our results indicate that one can use values of n of up to 10, 
without affecting significantly the calculated structural and dynamical magnitudes. 

PACS numbers: 71.15.Pd, 02.70.Ns, 31.15.Qg, 33.15.Vb 
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Molecular dynamics (MD) is a fundamental tool in 
atomistic materials simulation p|. A majority of prac- 
titioners have used classical, semiempirical interatomic 
potentials. This is necessary for the large sizes and long 
times required to simulate many processes of enormous 
scientific and technological interest, from materials de- 
formation and fracture 0] to protein folding A large 
effort has been devoted to develop interatomic potentials 
for many types of systems 4]. However, the quantita- 
tive reliability of such potentials in situations of bond 
formation and breaking is highly questionable. In such 
cases, it is imperative to use the much more expensive 
ab initio MD methods 0, @ , generally limited to a few 
hundred atoms and a few tens of picoseconds. Thus, it 
is essential to find methods that accelerate the integra- 
tion of the dynamical equations, thus allowing for longer 
simulations. 

In classical dynamics, one of such methods Q uses 
multiple time scales to integrate the equations of motion 
for systems with both fast and slow dynamical degrees 
of freedom. The same method can be used to compute 
separately the hard, short-ranged forces from the soft, 
long-ranged ones. De Vita and Car [8| have proposed to 
adapt 'on the fly' the parameters of a classical potential 
using sporadic or periodic evaluations of ab initio forces. 
In this work, drawing ideas of those previous works, we 
propose a simple method to speed up dramatically ab 
initio MD. In principle, it could be implemented by com- 
bining classical and ab initio forces. However, such an 
approach would still require to develop a suitable clas- 
sical force field for every new system with different in- 
teractions. Therefore, instead we take advantage of the 
fact that, while standard density functional forces require 
the simultaneous convergence of many parameters, much 
lower values of those parameters can still yield quite rea- 
sonable forces. Thus, by reducing drastically the size of 
the basis set, the Brillouin zone sampling, or the num- 



ber of selfconsistency iterations, it is possible to reduce 
the computer time by enormous factors, and still obtain 
forces which are considerably more reliable than those of 
classical interatomic potentials. 

To test our scheme, we have chosen the SIESTA 
method 0, 0, which is specially well suited to span 
the range from 'quick and dirty' calculations to fully con- 
verged ones. It uses density functional theory 

E3 (DFT), 

norm-conserving pseudopotentials 0] and a basis set of 
numerical atomic orbitals of strictly finite range [T3. fl4| . 
To calculate converged forces, we might typically use 
the generalized gradient approximation (GGA) to ex- 
change and correlation (with spin polarization if re- 
quired), double-^ polarized (DZP) basis orbitals with a 
relatively long range, fine integration grids in real and re- 
ciprocal space, and a well converged selfconsistency be- 
tween density and potential. For the cheap forces we 
may save on many different parameters, depending on 
the system and the properties studied. Thus, we may 
use the local density approximation (LDA), a minimal 
single-^ basis set with short range, a coarser integration 
grid in real space, just the T point in reciprocal space, 
and the non-selfconsistent Harris functional All to- 
gether, the cheap forces are typically hundreds of times 
faster to compute than the converged ones, and therefore 
they add a neglegible cost to the overall calculation, thus 
making unneccessary to resort to classical force fields. 

As usual, we use the Born-Oppcnheimcr approxima- 
tion and we treat the nuclei as classical particles, sub- 
ject to the Hellmann-Feynman forces (including all Pu- 
lay corrections). The equations of motion are solved with 
the standard velocity- Verlet algorithm what ensures 
the time reversibility of the trajectories |jfl. The atomic 
forces at time t are defined as F f ast (t) + AF(i), where 



AP(t) 



n(F conv (t) - F fast (t)) if (t/At mod n) = 
otherwise. 

(1) 
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Thus, the expensive converged forces F conv need to be 
evaluated only once every n time steps At. In those 'cor- 
rection steps', the trajectories generated by the cheap 
(fast) forces F f ast are corrected by applying a force 'kick' 
equal to the difference between the converged and fast 
forces at that time, multiplied by n. The factor n ac- 
counts for the concentration of the continuous force cor- 
rection F conv (t) — F f „, „t (t) in one out of every n steps. 
The method of Ref.lllJ based on the position- Ver let algo- 
rithm, was reported to have a better numerical stability 
in response to the correction 'kicks'. The efficiency of 
that method in the present context will be studied in 
future works. 

Figure n shows schematically the positions, velocities 
and forces of a particle moving in one dimension, gener- 
ated with our mixed-force algorithm. For simplicity, we 
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FIG. 1: Schematic a) position, b) velocity, and c) force, in 
arbitrary units, for a particle moving in one dimension, gen- 
erated by the mixed-force algorithm with zero initial velocity 
and a force correction interval of 10 time steps. For simplicity, 
the converged force is equal to zero and the cheap (fast) force 
is a negative constant. The periodic force correction 'kicks' 
are positive and invert the velocity. Notice that the position 
and velocity at the correction steps are equal to their correct 
converged value (zero). 

take the converged force and the initial velocity equal to 
zero, so that the correct converged position is also zero 
at all times. The mixed-force trajectory, for a constant 
negative fast force, shows periodic force kicks that change 
discontinuously the velocity and invert the trajectory at 
the correction steps. 

We have applied this method to simulate a system of 
64 silicon atoms at an average temperature of ~ 2000 K 
and an average pressure close to zero. This high temper- 
ature was intentionally chosen to test the method under 
specially stringent conditions, whith high kinetic energies 
and frequent formation and breaking of bonds. The sim- 
ulations were performed with the SIESTA program [loj| 
but standard Hamiltonian diagonalizations were used in- 
stead of order- TV methods 0,0] , because of the metallic 
character of liquid silicon. For the cheap forces we use 



the Harris functional, a minimal basis set with a range of 
3.5 and 4.0 Bohr for s and p orbitals, a real-space inte- 
gration grid with a plane wave cutoff of 40 Ry, and only 
the r fc-point. For the converged forces, we use the self- 
consistent Kohn-Sham functional in the LDA, a double-^ 
polarized (DZP) basis set with a range of 5.4, 6.5, and 
3.8 for s,p and d basis orbitals, a real-space grid with a 
80 Ry plane wave cutoff, and only the T fc-point. The 
forces are corrected according to equation every n 
time steps, with At = 1 fs. 

Figure [3 compares the magnitudes of the fast and con- 
verged forces, and of their difference. It can be seen that 
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FIG. 2: Decomposition of the total converged forces F conv 
into a cheaply evaluated component F/ ast , and a remainder 
F con „ — F/ ast . Represented are the average norms as a func- 
tion of time. The trajectory was generated for 64 Si atoms at 
2000 K, using the converged forces. 

the latter is a relativelly small and smooth correction, 
what explains why it may be evaluated and applied less 
frequently. 

Figure |31 represents the divergence of the trajectories, 
generated with different values of n, away from the ref- 
erence converged trajectory (which corresponds to n = 1 
in Eq. QJ). It can be seen that even the n — 10 tra- 
jectory diverges much more slowly than that obtained 
purely from the fast forces (labelled n = oo, i.e. with no 
force corrections). 

Figure 01 shows the total energy as a function of time. 
The energy conservation is considerably worse in the self- 
consistent converged-force trajectory (n = 1) than in the 
Harris-force trajectory (n — oo). This probably reflects 
larger effects of charge sloshings and analytic disconti- 
nuities in the forces, due to frequent level crossings in 
this highly disordered system. However, it is important 
to notice that the energy conservation in the mixed-force 
trajectory (n — 10) is similar to that in the converged 
trajectory. 

Despite the high simulation temperature, the charge 
transfer in elemental liquid silicon may be expected to 
be considerably smaller than in an ionic system, mak- 
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FIG. 3: Divergence of the mixed- force MD trajectories from 
the converged-force trajectory for the liquid silicon system. 
Ax is the average distance in atomic positions between the 
given and the reference trajectories. Force corrections were 
made every n — 2, 5, 10, 20, and oo (only fast forces) time 
steps. 
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FIG. 4: Total energy per atom, as a function of time, for 
the liquid silicon system. In the mixed-force (n = 10) and 
converged-force (n = 1) trajectories, the total energy was 
calculated at the correction steps, as the sum of the Kohn- 
Sham energy plus the nuclear repulsion and kinetic energies. 
In the fast-force trajectory (n = oo), it was calculated at 
every step, using the Harris-functional for the electronic part. 
The standard deviations are 0.9, 1.3, and 0.3 meV/atom for 
the n = 1, 10, and oo trajectories, respectively. 



ing the non-selfconsistent Harris functional specially ad- 
equate. In fact, we have seen that structural magnitudes 
like the bond length and bond angle distributions are not 
very different using the Kohn-Sham and Harris function- 
als. Therefore, we have also studied a more challenging 
system, liquid silica, using 72 atoms at a high average 
temperature of 5500 K and a low density of 0.42 g/cm 3 , 
typical of porous silica aerogels The distributions of 
bond legths and angles, presented in Figures |S] and arc 



indeed very different using the two functionals (n = 1 
and n — oo). Despite this, the mixed- force method, 
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FIG. 5: Radial pair distribution functions of liquid Si02 at 
5500 K, using the new method with different values of n. 
The converged forces, used for n = 1 and for the correction 
steps of n — 5 and 10, were obtained with a double £ plus 
polarization basis set, a real-space grid with a plane wave 
cutoff of 200 Ryd and only the F fc-point. The fast forces (that 
yield the n — oo curves when uncorrected) were calculated 
with a minimal basis set (single £), a real-space grid with a 
plane wave cutoff of 150 Ryd and only the F fc-point. 

with up to n — 10, yields essentially the same distribu- 
tions as the converged Kohn-Sham trajectory. Similar 
results, to be presented elsewhere, were obtained for an 
even more ionic system, liquid magnesium oxide, with 54 
atoms at 6500 K and 30 GPa. 

It might be expected that dynamical magnitudes are 
more sensitive than thermodynamic averages to changes 
in how the MD trajectories are obtained. Figure [7| shows 
the velocity autocorrelation function for the three sys- 
tems studied, as a function of the interval n between 
force corrections. As expected, the trajectory of the 
non-selfconsistcnt Harris functional is reasonably accu- 
rate only for elemental liquid silicon. But, in every case, 
the mixed-force method, with up to n — 10, yields es- 
sentially the same velocity autocorrelations as the con- 
verged Kohn-Sham trajectories. We have also calculated 
sclf-difussion coefficients from the average quadratic dis- 
tances traversed as a funtion of time. Thus, for liquid 
silicon we obtain, respectively, (2.4 ± 0.1), (2.5 ± 0.1), 
(2.6 ± 0.1), (2.6 ± 0.1), and (2.0 ± 0.1) x 10~ 4 cm 2 /s, 
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FIG. 6: Bond-angle distribution functions of the liquid SiCh 
system, using the new method with different values of n. 
The bond cutoff radii were chosen 10% larger than the first 
maximum in the corresponding pair distribution functions of 
Fig.0 
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FIG. 7: Velocity autocorrelation function Z(t) — 
(vi(t')vi(t' + i))/(vf (t')), as a function of the interval n be- 
tween force corrections, for a) liquid Si at 2000 K and zero 
pressure, b) liquid MgO at 6500 K and 30 GPa, and c) liquid 
Si0 2 at 5500 K and a density of 0.42 g/cm 3 . 

for n = 1,5, 10, 20 and oo. Again, the mixed-force value, 
even with n = 20, is the same, within the statistical er- 
ror, as that of the converged trajectory. This indicates 



that dynamical and kinetic magnitudes, as well as struc- 
tural or thermodynamic averages, are well reproduced 
even with quite large values of the boost factor n. 



In conclusion, we have presented a new method to 
greatly accelerate ab initio molecular dynamics simula- 
tions by combining cheap force evaluations with accurate 
converged ones. Our results show that the method is very 
robust with respect to the reduced accuracy of the cheap 
forces. Although the acceleration factor will undoubtly 
depend on the system simulated, our present results in- 
dicate that factors of 10 can be expected in most cases. 
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